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3D simulations of wind-jet interaction in massive X-ray binaries 
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ABSTRACT 



Context. High-mass microquasars may produce jets that will strongly interact with surrounding stellar winds on binary 
system spatial scales. 

Aims. We study the dynamics of the collision between a mildly relativistic hydrodynamical jet of supersonic nature and 
the wind of an OB star. 

Methods. We performed numerical 3D simulations of jets that cross the stellar wind with the code Ratpenat. 
Results. The jet head generates a strong shock in the wind, and strong recollimation shocks occur due to the initial 
overpressure of the jet with its environment. These shocks can accelerate particles up to TeV energies and produce 
gamma-rays. The recollimation shock also strengthens jet asymmetric Kelvin-Helmholtz instabilities produced in the 
wind/jet contact discontinuity. This can lead to jet disruption even for jet powers of several times 10 36 erg s _1 . 
Conclusions. High-mass microquasar jets likely suffer a strong recollimation shock that can be a site of particle acceler- 
ation up to very high energies, but also eventually lead to the disruption of the jet. 
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1. Introduction 

Jets of X-ray binaries (microquasars) are produced close to 
the compact object (black hole or neutron star) via ejection 
of material accreted from the stellar companion. Jet syn- 
chrotron emission was extensively observed (Ribo 2005), 
and its phenomenological properties and connections with 
other energy bands were thouroughly analyzed (Fender et 
al. 2004). The occurrence of collisionless shocks in micro- 
quasar jets can lead to efficient particle acceleration (Rieger 
et al. 2007) and non-thermal emission of synchrotron and 
inverse Compton origin and, possibly, from proton-proton 
collisions (see Bosch-Ramon & Khangulyan 2009 and ref- 
erences therein). To understand the dynamics of these jets, 
numerical simulations of their propagation were performed 
(Peter & Eichler 1995; Velazquez & Raga 2000; Perucho & 
Bosch- Ramon 2008; Bordas et al. 2009) together with sev- 
eral analytical treatments (Heinz & Sunyaev 2002; Heinz 
2002; Kaiser et al. 2004; Nalewajko & Sikora 2009; Araudo 
et al. 2009). Most of these works were valid for microquasars 
in general, although jets of high-mass systems require a spe- 
cific treatment due to the strong stellar wind. 

Perucho & Bosch-Ramon (2008) -PBR08 from now on- 
studied how the strong wind of an OB star can influence 
the jet dynamics at scales similar to the orbital separa- 
tion (~ 0.2 AU). Simulations in two dimensions of a hy- 
drodynamical jet interacting with an homogeneous (i.e. not 
clumpy) stellar wind were performed in cylindrical and slab 
symmetries. The results showed that a strong recollimation 
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shock is likely to occur at jet heights ~ 10 12 cm, which could 
lead to efficient particle acceleration and gamma-ray emis- 
sion, which was detected in several high-mass X-ray binaries 
(HMXB) (LS 5039, LS I +61 303, Cygnus X-l, Cygnus X- 
3; Aharonian et al. 2005, Albert et al. 2006, Albert et al. 
2007, Tavani et al. 2009). It was also found that jet dis- 
ruption could occur for jet kinetic luminosities as high as 
Lj ~ 10 36 erg s _1 because of jet instabilities produced by 
the strong and asymmetric wind impact. Such a Lj -value 
is ~ 0.1 — 1% of the Eddington luminosity, which is typi- 
cal for an X-ray binary with persistent jets (see Fig. 1 in 
Fender et al. 2003). All this means that the role of the stel- 
lar wind in high-mass microquasars (HMMQ) is not only 
to feed accretion, but it may be also relevant for particle 
acceleration within the binary system (i.e. in strong recolli- 
mation shocks) and a detectable jet at larger spatial scales. 
However, although the assumptions in PBR08 of a hydro- 
dynamical jet and a homogeneous wind seem reasonable 
at the involved jet heights of ~ 10 5 i?s c h (Sikora et al. 

2005) and under moderate clumping (see Owocki & Cohen 

2006) , the assumption of 2-dimensional symmetry is far less 
realistic. Therefore, 3-dimensional (3D) simulations of the 
wind impact on the jet are necessary to understand the 
non-thermal phenomena occurring in HMMQ. 

In this Letter, we present the results of 3D simulations 
of hydrodynamical jets interacting with a typical OB star 
wind. Wc find that even jets with Lj as high as several times 
10 36 erg s _1 may be disrupted, since the lower wind-jet 
momentum transfer due to wind sidewards escape is bal- 
anced by enhanced development of disruptive instabilities. 
The simulations also show a strong recollimation shock that 
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could efficiently accelerate particles. A deeper treatment of 
the radiative counterpart of the simulations presented here 
is in preparation (Bosch-Ramon, Khangulyan & Perucho, 
in preparation). 



2. Simulations 



We performed 3D simulations of two supersonic hydro- 
dynamical jets with Lj = 10 35 (Jet 1) and 10 37 erg s _1 
(Jet 2). The medium is an isotropic wind (as seen from 
the star) of a mass-loss rate of 10 -6 M yr _1 and a con- 
stant velocity of 2 x 10 8 cm s -1 (PBR08). The compan- 
ion star is assumed to be located at i? rb = 2 x 10 12 cm 
from the base of the jet in a direction perpendicular to 
the jet axis. Both jets are injected at a distance to the 
compact object of zq = 6 x 10 10 cm with an initial jet ra- 
dius Rj = 6 x 10 9 cm. Their initial velocity, temperature 
and Mach number are 0.55 c (1.7 x 10 10 cm s -1 ), 10 10 K 
and 17, respectively. The initial densities and pressures for 
the jets 1 and 2 are p\ — 0.088 p w and p2 = 8.8 p w , and 
Pj.i = 71 erg cm -3 and P h % = 7.1 x 10 3 erg cm -3 , respec- 
tively; p w — 3 x 10~ 15 g cm -3 is the wind density at zq. We 
used the coordinate z for the direction of the initial propa- 
gation of the jet, x for the direction connecting the jet base 
and the star, and y for the direction perpendicular to z and 
x. We assumed that the magnetic field has no dynamical 
influence in the evolution of the jet, the wind is continuous 
and homogeneous and the compact object is at the same 
orbital position during the simulation time, ~ 10 3 s, much 
smaller than the orbital period. 

We used a finite-difference code named Ratpenat, 
which solves the equations of relativistic hydrodynamics in 
three dimensions written in conservation form using high- 
rcsolution-shock-capturing methods. Ratpenat was paral- 
lelized with a hybrid scheme with both parallel processes 
(MPI) and parallel threads (OpenMP) inside each process 
(see Perucho et al. 2009). The simulations were performed 
in Mare Nostrum, at the Barcelona Supercomputing Centre 
(BSC) with up to 128 processors. The numerical grid box 
expands transversely 20 Rj on each side of the axis, which 
amounts to a total of 40 Rj, and 320 Rj along the axis. 
The numerical resolution of the simulation is of four cells 
per initial jet radius. This means that the final box is 
160 x 160 x 1280 cells. An extended grid is used in the 
transversal direction, composed by 80 cells with an increas- 
ing size, which brings the outer boundary 80 Rj farther from 
the axis. The resolution is low at Zg due to the amount of 
computational time needed to perform the simulations, but 
because jets initially expand, the effective resolution at the 
z of interest (~ 10 12 cm) in the main grid is ~ 16 cells per 
R 3 . To check possible resolution problems, a double reso- 
lution simulation of Jet 2 was performed up to one half of 
the jet final distance (which required 10 times longer). The 
accumulated error in this first portion of the simulation is 
negligible for our purposes in all the hydrodynamical vari- 
ables. 



3. Results 

Jet 1 propagates up to z < 1.6 x 10 12 cm after ~ 1.25 x 10 3 s. 
Just after injection, the difference between the velocity of 
the jet head and that of the backflow generates a low den- 
sity region around the jet base, and the backflow material 



fills it very slowly. The jet expands initially in this region, 
mixing with some of the cocoon material and finally gener- 
ating a thick shear layer with positive velocities. The back- 
flow interacts strongly with this outer region, generating 
instabilities that grow in the shear layer. These instabili- 
ties are quite asymmetric due to the nature of the cocoon 
in the plane of impact of the wind (see Fig.QJ. After expan- 
sion, the central region of the jet becomes underpressured 
with respect to its surroundings and recollimates until the 
formation of a reconfinement shock. Thism shock is quasi- 
steady, i.e. propagates very slowly from z s ~ 2 x 10 11 cm 
to 4 x 10 11 cm, as the pressure in the cocoon drops (see 
next section) . Downstream of the reconfinement shock, the 
instabilities, which were growing in the shear layer, propa- 
gate to the whole section of the jet as the internal jet flow 
is decelerated and becomes more sensitive to perturbations. 
This process ends up in the mixing and deceleration of the 
jet flow at z > 10 12 cm. 
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Fig. 1. Jet 1, cuts of axial velocity (in cms" 1 ) and rest- 
mass density (in gem -3 ) along the propagation axis. The 
upper two panels show a cut perpendicular to the star-jet 
plane, whereas the bottom panels show a cut in this plane. 
The scale of the figure respects the proportion with Fig. [31 



At the end of the simulation, the velocity of the bow 
shock at jet head is ~ 3 x 10 8 cm s . The fit for the bow- 
shock advance velocity results in a fast decrease with time: 
i>b s tx t~ 6 . Figure Q] shows two axial cuts of Jet 1 with 
axial velocity (two upper panels) and rest-mass density (two 
lower panels) , along the plane perpendicular to the star-jet 
plane (upper) and in this plane (lower) at the last snapshot. 
The deviation caused by the wind thrust can be observed in 
the plane xz, which is the plane formed by the optical star 
and the jet. The low velocity of the jet head at the end of 
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Fig. 2. Transversal cuts for the axial velocity, Mach number and tracer in Jet 1 at z ~ 1.3 x 10 12 cm. The Mach number 
is saturated for values higher than 20 (the wind). 



the simulation along with its final (destabilized) structure 
implies that the jet will be disrupted and will not propagate 
out of the binary system as a supersonic and collimated 
flow. Figure [2] shows transversal cuts of the axial velocity, 
Mach number and traceiQat z ~ 1.3x 10 12 cm. In this figure 
we can see the deformation of the bow shock caused by the 
wind thrust. The jet has been entrained up to its axis by the 
wind material (f max < 1). The maximum velocity in the 
jet fluid is still relatively fast (~ 1.5 x 10 10 cm s -1 ) despite 
the irregular morphology and mixing, but the average Mach 
number is close to one. 

The evolution of Jet 2 is similar to that of Jet 1 
from a qualitative point of view. It propagates up to z ~ 
2 x 10 12 cm in ~ 210 s. The jet expands more at the base 
because it is initially denser (more overpressured) , the ve- 
locity of the jet head is faster than in Jet 1, and the cocoon 
at this z is consequently filled more slowly. Therefore, the 
reconfincmcnt shock is stronger and occurs at larger z (see 
next section). The location of the shock changes with time 
from z ~ 6 x 10 11 cm to 10 12 cm (see Fig. [3]). The effect of 
the wind in the direction of the jet propagation is small, as 
seen in Figs. [3] and |H In the latter, the structure of the bow 
shock at z ~ 1.3 x 10 12 cm is observed to be more symmetric 
than that in Jet 1 (Fig. [2]). Figure [4] shows that the jet core 
is unmixed (/ = 1) and that the flow velocity is still as high 
as that in the injection point. At the end of the simulation, 
the velocity of the bow shock is ~ 8.4 x 10 9 cm s _1 , which 
is close to its initial speed (~ 9. x 10 9 cm s _1 ). However, 
the differences in the the physical conditions in the cocoon 
on both sides of the jet cause the development of asymmet- 
ric Kelvin-Helmholtz instabilities in the shear layer. Their 
propagation to the whole jet after the reconfinement shock 
triggers helical motions and distortions in the flow. 

At z <; Rorhi w e expect the density of the stellar wind to 
be oc z~ 2 . This implies a pressure drop and a jet adiabatic 
expansion. The helical instability triggered at z 5s i? rb 
will then develop in these changing conditions and may 
slow its growth (Hardee 2006). How much it will affect the 
stability of the jet flow in an expanding jet yet remains to 
be studied. However, since the instability reaches nonlinear 
amplitudes, the process is irreversible and we do not expect 
it to disappear. Therefore L } = 10 37 erg s^ 1 seems to be 
close to the minimum jet power to propagate out of the 
binary region without disruption. 
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1 The tracer, / = [0, 1], is a variable that indicates the compo- 
sition of the fluid in a given cell (see e.g., Perucho et al. 2005), 
with a value of corresponding to pure wind material, 1 to pure 
jet material, and any value between and 1 indicating the rela- 
tive amount of jet material in a cell where mixing has occurred. 
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Fig. 3. Same plots as in Fig. [TJ for Jet 2. 



4. Discussion 

The performed simulations show that a recollimation shock 
likely forms when the jet is crossing the wind. As men- 
tioned, this shock can trigger particle acceleration, and 
could also enhance instabilities that would destroy the jet. 
However, the jet- wind interaction simulated here is time- 
dependent, and the study of the stationary case calls for a 
specific simulation. Still, it is worthwhile to estimate here 
analytically whether the shock remains inside the binary 
region when the jet head is already outside. This will de- 
pend on whether the wind ram pressure can substitute the 
cocoon pressure to keep the recollimation shock inside the 
system. 

— 1/2 

From Falle (1991), we know that z s (x P c , with P c 
the pressure of the cocoon, if the jet power and mass flux 
are kept constant. Thus we need to know how this pres- 
sure changes with time. Following Scheck et al. (2002) and 
Perucho & Marti (2007), we find that P c oc t- 1 -"/ 2 for a ho- 
mogeneous ambient medium, where Vhs oc i Q ; and P c oc i~ 2 
when the ambient medium density decreases as p a oc z~ 2 . 
This implies z s oc i 1 / 2 + a / 4 and oc t, respectively. This result 
can be tested using our simulation, in which the ambient 
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Fig. 4. Same as in Fig. [5] for Jet 2 at z ~ 1.5 x 10 12 cm.The Mach number is saturated for values higher than 20 (the 
wind) . 



medium is still roughly homogeneous (at z Ss i?orb)- For 
Jet 2 the recollimation shock forms at z = 6 x 10 11 cm at 
t\ ~ 61 s. At t2 — 210 s, taking into account that the veloc- 
ity is basically constant (a = 0), z s (t2) / z s (t\) = (tn/ti) 1 ' 2 , 
which results in Zgfa) = 1.1 x 10 12 cm, in agreement with 
the simulation. We can then extrapolate to find out the 
time (£3) at which the shock would reach z = 2 x 10 12 cm. 
Taking t% — 210 s, and z s (t2) = 1-1 x 10 12 cm we obtain for 
p a oc z~ 2 : t 3 ~ 2i 2 ~ 400 s. 

We now calculate the time at which the pressure of the 
cocoon will fall below the pressure of the shocked wind. 
This is the time from which the recollimation shock will 
be determined by the latter, which is P w ~ 10 2 erg cm -3 
(see PBR08). Taking into account that the mean pressure 
in the cocoon is P c ~ 10 3 ergcm~ 3 at ti — 210 s, pres- 
sure equilibrium with the shocked wind will be reached at 
t ~ v/Pco/Pw t 2 3*2- This means that Jet 2 is around 
the limit to keep the recollimation shock inside the system 
by the wind ram pressure alone, under the simulated condi- 
tions. For hotter or denser jets, the recollimation shock can 
move out of the region of interest at a finite time, whereas 
for colder or lighter jets, this shock will stay inside the bi- 
nary system. This allows for continuous production of en- 
ergetic emission, but puts the jet in danger of disruption. 
An X-ray binary in which a disruption of the jet could be 
taking place is LS 5039 (Moldon et al. 2008), but higher 
angular resolution observations are required for a proper 
probe of the jet-wind interaction region. 

The scenario considered in this work could take place 
in several HMXB in the Galaxy. The luminosity function 
derived by Grimm et al. (2003) predicts ~ 3 HMXBs with 
L x = 10 35 erg/s. Following Fender et al. (2005), an HMXB 
with a 10 Mq black hole could produce a jet with a kinetic 
power between 10 35 and 10 38 erg/s, which is in the range 
of the simulations performed here. Although Grimm et al. 
(2003) do not offer a specific prediction for Lx < 10 35 erg/s, 
extrapolating the given luminosity function we deduce that 
there is room for ~ 10 sources in our Galaxy whose jets 
could be disrupted by the stellar wind. These objects may 
be bright at high energies, but faint or even quiet in radio. 

Several improvements are required to make 3D HMMQ 
jet simulations more realistic: to study a stationary case at 
the binary scales; to perform detailed calculations of the 
radiation produced in the jet-wind shocks; to introduce an 
inhomogeneous or clumpy stellar wind; to follow the jet 
head up to z > i? or b to study the effects of the decreasing 
wind density; and to account for the effect of the magnetic 
field on the dynamics of these jets, since the configuration 
of the magnetic field in the jet affects its stability (Hardee 
2007, Mizuno et al. 2007). All this is ongoing work to be 
presented elsewhere. 
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